In the appendices the forward method is compared with the backward
method using Neumann conditions.
Derivation of Neumann Boundary Conditions
Left Boundary
For the left boundary \(\frac{\delta
u}{\delta x} \Big|_{x=x_0} = 0\), let \(u_i = u(x_i,\tau)\) and use Taylor
expansion around \(x=x_0\)
\[\begin{align*}
u_0 &= u_0 \\
u_1 &= u_0
+ \Delta x \frac{\delta u}{\delta x}
+ \frac{1}{2} (\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ \frac{1}{6} (\Delta x)^3 \frac{\delta^3 u}{\delta x^3}
+ O\left(\left(\Delta x\right)^4\right) \\
u_2 &= u_0
+ 2\Delta x \frac{\delta u}{\delta x}
+ \frac{1}{2} (2\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ \frac{1}{6} (2\Delta x)^3 \frac{\delta^3 u}{\delta x^3}
+ O\left(\left(\Delta x\right)^4\right) \\
&= u_0
+ 2\Delta x \frac{\delta u}{\delta x}
+ \frac{4}{2} (\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ \frac{8}{6} (\Delta x)^3 \frac{\delta^3 u}{\delta x^3}
+ O\left(\left(\Delta x\right)^4\right) \\
&= u_0
+ 2\Delta x \frac{\delta u}{\delta x}
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ \frac{8}{6} (\Delta x)^3 \frac{\delta^3 u}{\delta x^3}
+ O\left(\left(\Delta x\right)^4\right)
\end{align*}\]
Eliminating the \(\frac{\delta^3 u}{\delta
x^3}\) terms
\[\begin{align*}
8u_1 - u_2 &= \left[8u_0
+ 8\Delta x \frac{\delta u}{\delta x}
+ \frac{8}{2} (\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ \frac{8}{6} (\Delta x)^3 \frac{\delta^3 u}{\delta x^3}
+ O\left(\left(\Delta x\right)^4\right)\right]
-\left[u_0
+ 2\Delta x \frac{\delta u}{\delta x}
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ \frac{8}{6} (\Delta x)^3 \frac{\delta^3 u}{\delta x^3}
+ O\left(\left(\Delta x\right)^4\right)\right] \\
&= 7u_0
+ 6\Delta x \frac{\delta u}{\delta x}
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ O\left(\left(\Delta x\right)^4\right)
\end{align*}\]
and eliminating the \(u_0\)
terms
\[\begin{align*}
8u_1 - u_2 -7u_0 &= \left[7u_0
+ 6\Delta x \frac{\delta u}{\delta x}
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ O\left(\left(\Delta x\right)^4\right)\right]
-7u_0 \\
&= 6\Delta x \frac{\delta u}{\delta x}
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ O\left(\left(\Delta x\right)^4\right)
\end{align*}\]
and knowing that each derivative is evaluated at \(x=x_0\) leaves
\[\begin{align*}
8u_1 - u_2 -7u_0 &= 6\Delta x \frac{\delta u}{\delta x}
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ O\left(\left(\Delta x\right)^4\right) \\
8u_1 - u_2 -7u_0 &= 6\Delta x (0)
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ O\left(\left(\Delta x\right)^4\right) \\
8u_1 - u_2 -7u_0 &= 2 (\Delta x)^2 \frac{\delta^2 u}{\delta x^2} +
O\left(\left(\Delta x\right)^4\right) \\
\frac{8u_1 - u_2 -7u_0}{2 (\Delta x)^2} &= \frac{\delta^2 u}{\delta
x^2} + O\left(\left(\Delta x\right)^4\right)
\end{align*}\]
From the heat equation
\[\begin{align*}
u_\tau - \kappa u_{xx} &= 0 \\
u_\tau &= \kappa u_{xx} \\
\frac{u_\tau}{\kappa} &= u_{xx}
\end{align*}\]
Applying this to the left boundary and letting \(\rho=\frac{\kappa \Delta \tau}{(\Delta
x)^2}\)
\[\begin{align*}
\frac{\delta^2 u}{\delta x^2} + O\left(\left(\Delta x\right)^4\right)
= \frac{8u_1 - u_2 -7u_0}{2 (\Delta x)^2}
&= \frac{u(x_0, \tau+\Delta \tau)-u(x_0, \tau)}{\kappa \Delta \tau}
= \frac{u_\tau}{\kappa} \\
\frac{1}{2}\frac{\kappa \Delta \tau}{(\Delta x)^2} (8u_1 - u_2 -7u_0)
&= u(x_0, \tau+\Delta \tau)-u_0 \\
\frac{1}{2} \rho (8u_1 - u_2 -7u_0) + u_0
&= u(x_0, \tau+\Delta \tau) \\
4 \rho u_1 - \frac{1}{2} \rho u_2 -\frac{7}{2}\rho u_0 + u_0
&= u(x_0, \tau+\Delta \tau) \\
\left(1-\frac{7}{2}\rho \right)u_0 + 4 \rho u_1 - \frac{1}{2} \rho u_2
&= u(x_0, \tau+\Delta \tau)
\end{align*}\]
This makes the elements of the first row of \(A_{\text{Explicit}}\)
\[\begin{align*}
A_{\text{Explicit}}[0,0] &= 1-\frac{7}{2}\rho \\
A_{\text{Explicit}}[0,1] &= 4 \rho \\
A_{\text{Explicit}}[0,2] &= - \frac{1}{2} \rho
\end{align*}\]
with the rest of the elements in the row equal to 0.
Right Boundary
For the right boundary \(\frac{\delta
u}{\delta x} \Big|_{x=x_N} = 0\), let \(u_i = u(x_i,\tau)\) and use Taylor
expansion around \(x=x_N\)
\[\begin{align*}
u_N &= u_N \\
u_{N-1} &= u_N
- \Delta x \frac{\delta u}{\delta x}
+ \frac{1}{2} (\Delta x)^2 \frac{\delta^2 u}{\delta^2 x}
- \frac{1}{6} (\Delta x)^3 \frac{\delta^3 u}{\delta^3 x}
+ O\left(\left(\Delta x\right)^4\right) \\
u_{N-2} &= u_N
- 2\Delta x \frac{\delta u}{\delta x}
+ \frac{1}{2} (2\Delta x)^2 \frac{\delta^2 u}{\delta^2 x}
- \frac{1}{6} (2\Delta x)^3 \frac{\delta^3 u}{\delta^3 x}
+ O\left(\left(\Delta x\right)^4\right) \\
&= u_N
- 2\Delta x \frac{\delta u}{\delta x}
+ \frac{4}{2} (\Delta x)^2 \frac{\delta^2 u}{\delta^2 x}
- \frac{8}{6} (\Delta x)^3 \frac{\delta^3 u}{\delta^3 x}
+ O\left(\left(\Delta x\right)^4\right) \\
&= u_N
- 2\Delta x \frac{\delta u}{\delta x}
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta^2 x}
- \frac{8}{6} (\Delta x)^3 \frac{\delta^3 u}{\delta^3 x}
+ O\left(\left(\Delta x\right)^4\right)
\end{align*}\]
Eliminating the \(\frac{\delta^3
u}{\delta^3 x}\) terms
\[\begin{align*}
8u_{N-1} - u_{N-2} &= \left[8u_N
- 8\Delta x \frac{\delta u}{\delta x}
+ \frac{8}{2} (\Delta x)^2 \frac{\delta^2 u}{\delta^2 x}
- \frac{8}{6} (\Delta x)^3 \frac{\delta^3 u}{\delta^3 x}
+ O\left(\left(\Delta x\right)^4\right)\right]
-\left[u_0
- 2\Delta x \frac{\delta u}{\delta x}
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta^2 x}
- \frac{8}{6} (\Delta x)^3 \frac{\delta^3 u}{\delta^3 x}
+ O\left(\left(\Delta x\right)^4\right)\right] \\
&= 7u_N
- 6\Delta x \frac{\delta u}{\delta x}
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta^2 x}
+ O\left(\left(\Delta x\right)^4\right)
\end{align*}\]
and eliminating the \(u_0\)
terms
\[\begin{align*}
8u_{N-1} - u_{N-2} -7u_N &= \left[7u_N
- 6\Delta x \frac{\delta u}{\delta x}
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta^2 x}
+ O\left(\left(\Delta x\right)^4\right)\right]
-7u_N \\
&= 6\Delta x \frac{\delta u}{\delta x}
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta^2 x}
+ O\left(\left(\Delta x\right)^4\right)
\end{align*}\]
and knowing that each derivative is evaluated at \(x=x_N\) leaves
\[\begin{align*}
8u_{N-1} - u_{N-2} -7u_N &= 6\Delta x \frac{\delta u}{\delta x}
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta^2 x}
+ O\left(\left(\Delta x\right)^4\right) \\
8u_{N-1} - u_{N-2} -7u_N&= 6\Delta x (0)
+ 2 (\Delta x)^2 \frac{\delta^2 u}{\delta^2 x}
+ O\left(\left(\Delta x\right)^4\right) \\
8u_{N-1} - u_{N-2} -7u_N &= 2 (\Delta x)^2 \frac{\delta^2
u}{\delta^2 x} + O\left(\left(\Delta x\right)^4\right) \\
\frac{8u_{N-1} - u_{N-2} -7u_N}{2 (\Delta x)^2} &= \frac{\delta^2
u}{\delta^2 x} + O\left(\left(\Delta x\right)^4\right)
\end{align*}\]
From the heat equation
\[\begin{align*}
u_\tau - \kappa u_{xx} &= 0 \\
u_\tau &= \kappa u_{xx} \\
\frac{u_\tau}{\kappa} &= u_{xx}
\end{align*}\]
Applying this to the left boundary and letting \(\rho=\frac{\kappa \Delta \tau}{(\Delta
x)^2}\)
\[\begin{align*}
\frac{\delta^2 u}{\delta^2 x} + O\left(\left(\Delta x\right)^4\right)
= \frac{8u_{N-1} - u_{N-2} -7u_N}{2 (\Delta x)^2}
&= \frac{u(x_N, \tau+\Delta \tau)-u(x_N, \tau)}{\kappa \Delta \tau}
= \frac{u_\tau}{\kappa} \\
\frac{1}{2}\frac{\kappa \Delta \tau}{(\Delta x)^2} (8u_{N-1} - u_{N-2}
-7u_N)
&= u(x_N, \tau+\Delta \tau)-u_N \\
\frac{1}{2} \rho (8u_{N-1} - u_{N-2} -7u_N) + u_N
&= u(x_N, \tau+\Delta \tau) \\
4 \rho u_{N-1} - \frac{1}{2} \rho u_{N-2} -\frac{7}{2}\rho u_N + u_N
&= u(x_N, \tau+\Delta \tau) \\
- \frac{1}{2} \rho u_{N-2} + 4 \rho u_{N-1} + \left(1-\frac{7}{2}\rho
\right)u_N
&= u(x_N, \tau+\Delta \tau)
\end{align*}\]
This makes the elements of the first row of \(A_{\text{Explicit}}\)
\[\begin{align*}
A_{\text{Explicit}}[N,N-2] &= - \frac{1}{2} \rho \\
A_{\text{Explicit}}[N,N-1] &= 4 \rho \\
A_{\text{Explicit}}[N,N] &= 1-\frac{7}{2}\rho
\end{align*}\]
with the rest of the elements in the row equal to 0.
Between the Boundaries
Within the boundaries, utilize Taylor series about \(x=x_i\) for
\[\begin{align*}
u_{i-1} &= u_i
- \Delta x \frac{\delta u}{\delta x}
+ \frac{1}{2}(\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
- \frac{1}{6}(\Delta x)^3 \frac{\delta^3 u}{\delta x^3}
+ O\left(\left(\Delta x\right)^4\right) \\
u_i &= u_i \\
u_{i+1} &= u_i
+ \Delta x \frac{\delta u}{\delta x}
+ \frac{1}{2}(\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ \frac{1}{6}(\Delta x)^3 \frac{\delta^3 u}{\delta x^3}
+ O\left(\left(\Delta x\right)^4\right)
\end{align*}\]
Eliminating the \(\frac{\delta u}{\delta
x}\) and \(\frac{\delta^3 u}{\delta
x^3}\) terms
\[\begin{align*}
u_{i-1} + u_{i+1} &= \left[u_i
- \Delta x \frac{\delta u}{\delta x}
+ \frac{1}{2}(\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
- \frac{1}{6}(\Delta x)^3 \frac{\delta^3 u}{\delta x^3}
+ O\left(\left(\Delta x\right)^4\right)\right]
+ \left[u_i
+ \Delta x \frac{\delta u}{\delta x}
+ \frac{1}{2}(\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
+ \frac{1}{6}(\Delta x)^3 \frac{\delta^3 u}{\delta x^3}
+ O\left(\left(\Delta x\right)^4\right)\right] \\
&= 2u_i + (\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
\end{align*}\]
and eliminating the \(u_i\) term
leaves
\[\begin{align*}
u_{i-1} + u_{i+1} - 2u_i &= \left[2u_i
+ (\Delta x)^2 \frac{\delta^2 u}{\delta x^2}\right] - 2u_i \\
u_{i-1} + u_{i+1} - 2u_i&=(\Delta x)^2 \frac{\delta^2 u}{\delta x^2}
\\
\frac{u_{i-1} + u_{i+1} - 2u_i}{(\Delta x)^2} &= \frac{\delta^2
u}{\delta x^2}
\end{align*}\]
From the heat equation
\[\begin{align*}
u_\tau - \kappa u_{xx} &= 0 \\
u_\tau &= \kappa u_{xx} \\
\frac{u_\tau}{\kappa} &= u_{xx}
\end{align*}\]
Applying this and letting \(\rho=\frac{\kappa \Delta \tau}{(\Delta
x)^2}\) yields
\[\begin{align*}
\frac{\delta^2 u}{\delta^2 x} + O\left(\left(\Delta x\right)^4\right)
= \frac{u_{i-1} + u_{i+1} - 2u_i}{(\Delta x)^2}
&= \frac{u(x_i, \tau+\Delta \tau)-u(x_i, \tau)}{\kappa \Delta \tau}
= \frac{u_\tau}{\kappa} \\
\frac{\kappa \Delta \tau}{(\Delta x)^2} (u_{i-1} + u_{i+1} - 2u_i)
&= u(x_i, \tau+\Delta \tau)-u_i \\
\rho (u_{i-1} + u_{i+1} - 2u_i) + u_i
&= u(x_i, \tau+\Delta \tau) \\
\rho u_{i-1} + \rho u_{i+1} - 2 \rho u_i + u_i
&= u(x_i, \tau+\Delta \tau) \\
\rho u_{i-1} + (1 - 2 \rho) u_i + \rho u_{i+1}
&= u(x_i, \tau+\Delta \tau)
\end{align*}\]
This makes the elements of the rows of \(A_{\text{Explicit}}\) between the
boundaries
\[\begin{align*}
A_{\text{Explicit}}[i,i-1] &= \rho \\
A_{\text{Explicit}}[i,i] &= 1-2\rho \\
A_{\text{Explicit}}[i,i+1] &= \rho
\end{align*}\]
with the rest of the elements in each row equal to 0.
Convert from \(A_{\text{Explicit}}\) to \(A_{\text{Implicit}}\)
From the work in the previous sections, the \(A_{\text{Explicit}}\) is found to be
\[
A_{\text{Explicit}} =
\begin{bmatrix}
1-\frac{7}{2}\rho & 4p & -\frac{1}{2}\rho & 0 & 0 &
\cdots & 0 \\
\rho & 1-2\rho & \rho & 0 & 0 & \cdots & 0 \\
0 & \rho & 1-2\rho & \rho & 0 & \cdots & 0 \\
\vdots & & \ddots & \ddots & \ddots & & \vdots
\\
0 & \cdots & 0 & \rho & 1-2\rho & \rho & 0 \\
0 & \cdots & 0 & 0 & \rho & 1-2\rho & \rho \\
0 & \cdots & 0 & 0 & -\frac{1}{2}\rho & 4p &
1-\frac{7}{2}\rho
\end{bmatrix}
\]
Since the goal is to compare the use of Neumann versus Dirichlet
boundary conditions, and the method given for the Dirichlet boundary
conditions was implicit, \(A_{\text{Implicit}}\) was found using \(A_{\text{Implicit}}=2I-A_{\text{Explicit}}\)
\[\begin{align*}
A_{\text{Implicit}}&=2I-A_{\text{Explicit}} \\
&= 2
\begin{bmatrix}
1 & 0 & 0 & \cdots & 0 & 0 & 0 \\
0 & 1 & 0 & \cdots & 0 & 0 & 0 \\
0 & 0 & 1 & \cdots & 0 & 0 & 0 \\
\vdots & & & \ddots & & & \vdots \\
0 & 0 & 0 & \cdots & 1 & 0 & 0 \\
0 & 0 & 0 & \cdots & 0 & 1 & 0 \\
0 & 0 & 0 & \cdots & 0 & 0 & 1
\end{bmatrix}
-
\begin{bmatrix}
1-\frac{7}{2}\rho & 4p & -\frac{1}{2}\rho & 0 & 0 &
\cdots & 0 \\
\rho & 1-2\rho & \rho & 0 & 0 & \cdots & 0 \\
0 & \rho & 1-2\rho & \rho & 0 & \cdots & 0 \\
\vdots & & \ddots & \ddots & \ddots & & \vdots
\\
0 & \cdots & 0 & \rho & 1-2\rho & \rho & 0 \\
0 & \cdots & 0 & 0 & \rho & 1-2\rho & \rho \\
0 & \cdots & 0 & 0 & -\frac{1}{2}\rho & 4p &
1-\frac{7}{2}\rho
\end{bmatrix} \\
&=
\begin{bmatrix}
2 & 0 & 0 & \cdots & 0 & 0 & 0 \\
0 & 2 & 0 & \cdots & 0 & 0 & 0 \\
0 & 0 & 2 & \cdots & 0 & 0 & 0 \\
\vdots & & & \ddots & & & \vdots \\
0 & 0 & 0 & \cdots & 2 & 0 & 0 \\
0 & 0 & 0 & \cdots & 0 & 2 & 0 \\
0 & 0 & 0 & \cdots & 0 & 0 & 2
\end{bmatrix}
-
\begin{bmatrix}
1-\frac{7}{2}\rho & 4p & -\frac{1}{2}\rho & 0 & 0 &
\cdots & 0 \\
\rho & 1-2\rho & \rho & 0 & 0 & \cdots & 0 \\
0 & \rho & 1-2\rho & \rho & 0 & \cdots & 0 \\
\vdots & & \ddots & \ddots & \ddots & & \vdots
\\
0 & \cdots & 0 & \rho & 1-2\rho & \rho & 0 \\
0 & \cdots & 0 & 0 & \rho & 1-2\rho & \rho \\
0 & \cdots & 0 & 0 & -\frac{1}{2}\rho & 4p &
1-\frac{7}{2}\rho
\end{bmatrix} \\
&=
\begin{bmatrix}
2-\left(1-\frac{7}{2}\rho\right) & 0-4p &
0-\left(-\frac{1}{2}\rho\right) & 0 & 0 & \cdots & 0 \\
0-\rho & 2-\left(1-2\rho\right) & 0-\rho & 0 & 0 &
\cdots & 0 \\
0 & 0-\rho & 2-\left(1-2\rho\right) & 0-\rho & 0 &
\cdots & 0 \\
\vdots & & \ddots & \ddots & \ddots & & \vdots
\\
0 & \cdots & 0 & 0-\rho & 2-\left(1-2\rho\right) &
0-\rho & 0 \\
0 & \cdots & 0 & 0 & 0-\rho & 2-\left(1-2\rho\right)
& 0-\rho \\
0 & \cdots & 0 & 0 & 0-\left(-\frac{1}{2}\rho\right)
& 0-4p & 2-\left(1-\frac{7}{2}\rho\right)
\end{bmatrix} \\
A_{\text{Implicit}}
&=
\begin{bmatrix}
1+\frac{7}{2}\rho & -4p & \frac{1}{2}\rho & 0 & 0 &
\cdots & 0 \\
-\rho & 1+2\rho & -\rho & 0 & 0 & \cdots & 0 \\
0 & -\rho & 1+2\rho & -\rho & 0 & \cdots & 0 \\
\vdots & & \ddots & \ddots & \ddots & & \vdots
\\
0 & \cdots & 0 & -\rho & 1+2\rho & -\rho & 0 \\
0 & \cdots & 0 & 0 & -\rho & 1+2\rho & -\rho \\
0 & \cdots & 0 & 0 & \frac{1}{2}\rho & -4p &
1+\frac{7}{2}\rho
\end{bmatrix}
\end{align*}\]
Backward (Implicit) Method with Neumann Boundary Conditions
Code to set up the \(A_{\text{Implicit}}\) Matrix
# Establish the A matrix for implicit scheme with Neumann boundary conditions
def setUpNeumannA(N, rho):
# Initialize a (N+1) by (N+1) matrix of 0s for A
A = np.zeros((N+1, N+1))
# Generate the first row of matrix A from the boundary conditions
A[0][0:3] = [1+((7/2)*rho), (-4*rho), (0.5*rho)]
# Generate the last row of matrix A from the boundary conditions
A[N][N-2:N+1] = [(0.5*rho), (-4*rho), (1+(7/2)*rho)]
# Fill in the rest of the A matrix as expected for the implicit method
for i in range(0,N-1):
j = i+1
A[j][j-1:j+2] = [-rho, (1+(2*rho)), -rho]
return(A)
Main Neumann Function for Backward (Implicit) Method
def Neumann(Stock_or_Rate, M, plot=True):
# M = number of simulations/iterations - if M is large enough it will go very little ahead in time for each iteration
x, u = Stock(Stock_or_Rate) # Stock Prices
#x, u = Rate(Stock_or_Rate) # Exchange Rates
res = np.zeros((len(u), M)) # Initialize matrix res, which holds the vector of stock prices generated for each iteration
# Set N equal to 1 less than the length of the u-vector
N = len(u)-1
dX = (b-a)/N # the distance between each closing price - i.e. time between daily closing prices
dT = T/M # the distance in time between each iteration
rho = kappa*dT/(dX*dX)
A = setUpNeumannA(N, rho)
if plot:
plt.figure(figsize=(20,10))
plt.plot(x, u)
for j in range(M):
u = np.linalg.solve(A, u)
res[:,j] = u # Store the vector of stock prices generated for this iteration in column j of matrix res
if plot:
plt.plot(x, u)
plt.grid(True)
plt.xlabel('$x$')
plt.ylabel('$u(x,\tau)$')
if plot:
plt.show()
return(res[:,M-1]) # Return the last column of matrix res (The vector of prices from the final iteration)
Effect of Changing Kappa \(\left( \kappa
\right)\) for Backward (Implicit) Method with Neumann Boundary
Conditions
As with Dirichlet conditions, increasing the value of \(\kappa\) enhances smoothing; however, using
too large a \(\kappa\) can cause too
much smoothing to the point that the heat equation smoothing converges
to a straight, horizontal line corresponding to the average price of the
price curve. On the other hand, using too small a value for \(\kappa\) will lead to very little
smoothing, to the point that the smoothed curve will converge to the
true stock price curve.
The following plots show the heat equation smoothing of the Credit
Suisse stock price using \(\kappa\)
values of 0.00001, 0.0001, 0.001, 0.01, 0.1, 1, 10, 100, and 1000,
holding all other parameters constant and using only 100 iterations for
each.
\(\kappa\) = 0.00001
kappa = 0.00001
ktn9 = Neumann(CS, 100)

\(\kappa\) = 0.0001
kappa = 0.0001
ktn8 = Neumann(CS, 100)

\(\kappa\) = 0.001
kappa = 0.001
ktn7 = Neumann(CS, 100)

\(\kappa\) = 0.01
kappa = 0.01
ktn6 = Neumann(CS, 100)

\(\kappa\) = 0.1
kappa = 0.1
ktn5 = Neumann(CS, 100)

\(\kappa\) = 1
kappa = 1
ktn4 = Neumann(CS, 100)

\(\kappa\) = 10
kappa = 10
ktn3 = Neumann(CS, 100)

\(\kappa\) = 100
kappa = 100
ktn2 = Neumann(CS, 100)

\(\kappa\) = 1000
kappa = 1000
ktn1 = Neumann(CS, 100)
